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ABSTRACT 

We analyze damping of oscillations of general relativistic superfluid neutron stars. 
To this aim we extend the method of decoupling of superfluid and normal oscilla- 
tion modes first suggested in [Gusakov & Kantor PRD 83, 081304(R) (2011)]. All 
calculations are made self-consistently within the finite temperature superfluid hy- 
drodynamics. The general analytic formulas are derived for damping times due to 
the shear and bulk viscosities. These formulas describe both normal and superfluid 
neutron stars and are valid for oscillation modes of arbitrary multipolarity. We show 
that: (i) use of the ordinary one-fluid hydrodynamics is a good approximation, for 
most of the stellar temperatures, if one is interested in calculation of the damping 
times of normal /-modes; (ii) for radial and p-modes such an approximation is poor; 
(Hi) the temperature dependence of damping times undergoes a set of rapid changes 
associated with resonance coupling of neighboring oscillation modes. The latter effect 
can substantially accelerate viscous damping of normal modes in certain stages of 
neutron-star thermal evolution. 
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1 INTRODUCTION 



circumferential radius R 

-i!4 „ „™-3 



10 km, and the central 



Neutron stars (NS) are compact objects with the mass M ~ 

density p c several times higher than the nuclear density po ~ 2.8 x 10 14 g cm" 3 . They are interesting because of extreme 
conditions in their interiors and a wide variety of associated astrophysical phenomena. In particular, internal instabilities or 
external perturbations ca n excite NS oscillations, whic h are potentially detectable by the next-generation gravitational wave 



interferometers (see, e.g. 



i excite ins oscillations, wmcn arc potentially detectable by tn 
Andersson fc KokkotaJboOll : lAnderssonlbood : lOwerJiioim . It 



is very probable, that quasiperiodic 



oscillations of ele ctromagnetic radiation observed in the tails of the giant gamma-ray flares are connected with oscillations 



in NS crust (e.g.. Ilsrael et alj|2005 l: Istrohmaver fc Wattsll2005l , bood : IWatts fc StrohmaverlbOQTj). and that seismology would 



V o , . .. . , .. . . 1 I 

becom e a significant source of information about NSs in the nearest future (jAbbot et al.l 
201 lh 



2007; Watts 2011; Andersson et al 



For the correct interpretation of already existing and future observations one requires a well-developed theory of oscillating 
NSs. It should, in particular: (i) be based on the general relativity theory, since NSs are relativistic objects; (ii) employ an 
adequate model of superdense matter, including realistic equation of state and parameters of baryon superfluidity; (Hi) 
correctly account for the effects of baryon superfluidity on the hydrodynamics of NS matter. 



Let us discuss briefly a (key) role of superfluidity. According to numerous microscopic calculations (see, e.g.. lLombardo fc Schulze 



200 



3), baryon matter in the internal layers of neut ron stars becomes superfl uid at T < 10 8 -10 10 K. It is very difficult to interpret 



the observational data on pulsar glitches (see, e.g. , Chamel fc Haensell2008l ) and cooling of NSs ( Yakovlev. Levenfish fc Shibanov 
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1999l ; lYakovlev fc Pethici3l20o4) without invok ing baryon superfluidity. Recent real-time observations of cooling NS in Cassio- 
pea A supernova remnant ( Heinke fc Ho 2010h also present a strong argum ent in favor of the exi stence of baryon superfluidity 
in the NS core. The observations w e re exp laine d bvlShternin et al.l ( 2011! ) and lPage et alj ( 201 In within a scenario, suggested 



for the first time in Gu sakov et al . (200 



t cn max 



and lPage et al.l (J200J), and assuming mild neutron superfluidity (with maximum 
7 -f- 9 x 10 8 K) and strong proton superconductivity (with maximum proton critical 



neutron critical temperatures T ( 
temperatures T cp max >2-f3x 10 9 K) in the NS core. 

Combined analysis of all the three factors (i)-(m) is a formidable task for the oscillation theory so in the literature 



Chandrasekharl (11964) and Thorne & Camoolattaro (1967) and were 


further developed in many subsequent papers (see, e.g., 


loser & Thorne 1973: Detweiler & iDserl 


1973: Lindblom & Detweiler 198 


i Detweiler & Lindblom 19851: Cutler & Lindbloml 


1987; 


Cutler. Lindblom & Splinter! 199C 


: Chandrasekhar & Ferrari 


1991 


: Kokkotas & Schutzl 199a" Yoshida & Leel l2003al: 


Lin. Andersson & ComerlboOSl and a review oflKokkotas & Schmidt Il999^ 


). When studying oscillations of NSs, most of these 



works considered ordinary one-fluid relativistic hydrodynamics. 

Meanwhile it is well known that superfluidity leads to appearance of additional velocity fields , desc ribing the superfluid 



degrees of freedom (e.g., Khalatnikov fc Lebedev 19821 : Khalatnikov 



complicates the hydrodynamics of NS matter, making it multi-fluid i|Mendel]||l991al lbl; iGusakov fc Andersson! bOQcJ ). In ad 



1989 



Carter fc Khalatnikov 1992). This subs tantially 



dition, superfluidity affects the kinetic coefficients (such as bulk and she ar viscosities) and also requires additional viscous 
coefficients to be introduced (see Gusakov 2007 ; Gusakov fc Kantor 2008! for details) . 



Oscillations of superfluid NSs have been studied actively only in the last two decades (see, e.g..lLeell995l:lLindblom fc Mendel] 



2000 



2009 



Prix & Ricutord 2002; Yoshida & Lcc 2003b; Prix, Comer & Andersson 2004; Samuelsson & Andersson 2009; Wong, Lin & Leung 



Passamonti fc Anderssonll2011 



2012), starting from the pioneering papers by lEpsteinl 1 19881) and lLindblom fc Mendelll 
( 1994]). However, most of these works neglect general relativity effects and employ zero temperature (T = 0) limit of superfluid 
hydrod ynamics (i.e., hydrodynamics, applicable only at T = 0). Within the general relativity theory os c illations were discussed 



only b y Comer 



(2006) 



Lin et al 



Andersson et al 



janglqis fc Lin Jl999h: Andersson. Comer fc Langlois J2OO2I): Yoshida fc Lee! ( 2003a ) ; IGusakov fc Andersson! 
(|2008!) ; lKantor fc Gusakovf (|201lh : Chugunqv fc Gusakov! d201ll). but most of thes e works used zero tempera- 



2002; Lin et al"1l2008l ) the presence of superfluid 



ture hydrodynamics. Moreover, in some of these papers (e.g., 
component was modeled by an artificial (polytropic) equation of state which does not represent any specific microphysical 
model. 



Only in the recent papers Gusakov fc Andersson! ( 20061 ) : Kantor fc Gusakov I 201ll ); Chugunov fc Gusakov ( 201 lh an 
attempt was made to self-consistently calculate the oscillation spectra using a realistic model of superdense matter and 
allowing for the effects of finite stellar temperatures, ft was shown that in many cases an approximation T = is n ot justified 
and, moreover, it can lead to qualitatively incorrect results (jKantor fc Gusakovll201ll ; IChugunov fc G usakov 2 01 lh ■ 



Of particular interest is the question of how superfluidity influences dissipation of neutron star oscillations. It is of 
extreme importance, for instance, for understanding physical conditions under which a rotating NS becomes unstable with 
respect to excitation of vario us oscillations (e.g., r-modes), and for estimating gravitational radiation from such stars (e.g., 
Andersson fc Kokkotasll200ll ). 



There were several serious and successful attempts to allow for the effects of superfluidity when studying the dissipation 
of oscillations in NSs (see, e.g ., Lindblom fc Mendelll! 1 995. 200 0; iLee fc Y oshida 2003; Haskell, Ander sson fc Passam onti 200S 



, ; _ 7 , ^ — ■— ^^ji ■ p^^^^^ _ = ^ii 1 

Andersson. Glampedakis fcHaskel l 2009; Haskell fc Anderssonll2010l ; IPassamonti fc Glampedakisll2012l ). but all 0} them con 



sidered Newtonian stars and used the T = superfluid hydrodynamics. The self -consistent analysis of dis sipation in superfluid 
NSs was only recently performed for a simple case of a radially oscillating NS (jKantor fc Gusakovll201ll ). 

The aim of the present paper is to fill this gap and to consider, for the first time, dissipation of nonradial oscillations 
in general relativistic superfluid NSs employing realistic microphysics input with accurate treatment of the effects of finite 
stellar temperatures. 

The paper is organized as follows. Relativistic superfluid hydrodynamics is briefly reviewed in Sec. [2] Sec. [3] discusses an 
unperturbed star and introduces variables describing small deviations of NS from equilibrium. In Sec. [4] we derive expressions 
for the oscillation energy and its dissipation rates due to bulk and shear viscosities. In Sec. [5] the equations that govern 
oscillations of superfluid NSs are explicitly written out. Sec. [6] describes the approach to study dissipation of superfluid NS 
oscillations. This approach is applied for a detailed numerical analysis of realistic models of oscillating neutron stars in Sec. 
[7] Sec. [8] presents a summary of our results. 



In what follows, we use the system of units in which c = &b = 1, where c is the speed of light and fce is the Boltzmann 
constant. 
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2 DISSIPATIVE SUPERFLUID HYDRODYNAMICS 



In this paper we consider, for simplicity, npe-matter in NS cores, that is matter composed of neutrons (n), protons (p), and 
electrons (e). Because both protons and neutrons can be in the superfluid state, one has to use the relativistic hydrodynamics 
of superfluid mixtures to study oscillations of NSs. Here we briefly discuss the corresponding equation s to establish notations 
and t o make the prese ntation more se lf-contained. Our consider ation closely follows the papers by Gusakov fc Anderssonl 
(2006); iGusakovl 1(2007]) and, especially, iKantor fc Gusakov! (|2011f ) The reader is referred to these works for more details. 

The main distinctive feature of superfluid hydrodynamics is the presence of several velocity fields in the mixture. In 
our case, these are the four-velocity of the 'normal' (nonsuperfluid) component of matter (electrons and Bogoliubov 
excitations of neutrons and protons) as well as the 'four- velocities' of superfluid neutrons v^, a -. and superfluid protons v^ p y 
In what follows instead of the velocities v^, a -, and w^/pj it will be convenient to use the four- vectors = fii[v^u) — 
where \jh is the relativistic chemical potential for particle species i — n or p. A presence of several velocity fields modifies the 
expressions for the current densities of neutrons and protons jf p -. , 

+ Y ik w? h) (1) 



Hi 



in comparison with the standard expression — mu^. The electron current density has a standard form, 



J 



(<0 



(2) 



Here and below the subscripts i and k refer to nucleons: i, k = n, p; n; is the number density of particle species I = n, p, 
e. Unless otherwise stated the summation is assumed over the repeated nucleon indices i, k and over the spacetime indices 
p, v, ... (Greek letters). In Eq. ffl Yik is the relativistic entrainment matrix, which is a generalization of the concept of 
superfluid density (see, e.g., Khalatnikov 19891) to the c ase of relativistic mixtures. In the nonrelativistic theory, a similar 
matrix was first considered bv lAndreev fc Bashkinl (|l975T ). The matrix Yik is symmetric, Yik = Y^i, and is expressed in terms 



of the Landau parameters F* k of a symmetric nuclear matter and universal functions of temperature, <E>i, as described in 
Gusakov. Kantor fc Haensell (2009b). In beta-equilibrium it can be presented as a function of density p and the combinations 



T/T cn and T/T cp : Yik = Yik(p,T/T cn ,T/T cp ), where T is the temperature; T cn (p) and T cp (p) are the density-dependent 
neutron and proton critical temperatures, respectively. If, for example, T > T cn then all neutrons are normal. The important 
property of the matrix Yu, is that for any nonsuperfluid species I = n or p, the corresponding elements Yik of this matrix 
vanish. 

In the present paper we consider NS oscillations, whose frequencies are well below the electron and proton plasma 
frequencies. In that case the quasineutrality condition, n c — n p , should hold in an oscillating star, from which it follows (for 
a nonrotating non- magnetized NS) jf 1 > = j 1 ?, or, in view of flj and pjl. 



Y,_ 



P kW (k) 



(3) 



Below we assume that this condition is always satisfied. It relates the four- vectors and W( p y 

In what follows, along with u M and wf^ it will be convenient to introduce the quantity , describing superfluid degrees 
of freedom, as well as the quantity which we call the 'baryon four-velocity' U£s (notice, however, that it is not a four-velocity 



in the usual sense, because generally UK^U^ M 7^ — 1, see Eq. (|45p and the footnote SI below). They are defined by the formulas 



Y n kV) 



(k) 



U 



(bi 



Mb 



(4) 
(5) 



where rib = n n + n p is the baryon number density. Notice that, as follows from Eqs. {T}-((3]), the baryon current density 



■7(b) ~ J(n) + J(p) ^ S re l a ted to by the standard equation, 

3(b) ~ nh ^(b) ' 

while j?\ equals 



(6) 



(7) 



Together w ith the quasineutrality condition (n e = n p ) and Eq. §3$, the equations of superfluid hydrodynamics include 
( Gusakov! [20071 ): 

(i) Continuity equations for baryons (b) and electrons (e), 



- fi 



'(e); M 



= 0, 
= 0; 



(8) 
(9) 
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(ii) Energy-momentum conservation 
_ (p + g-j u m u " + + y ik (wfawfa + fit yj^ k) u v + fi k W(i)« M J 



U-y-S + US:-, 



9-yS u. t 



(Hi) Potentiality condition for superfluid motion of neutrons 

[W(n)n + (Mn + = df± [ltf(n> + (^n + X n )w„] , 

as well as (iv) the second law of thermodynamics 



de = T dS + y, e dn e + fn dru + d (wfy W^) a ) 



(10) 
(11) 
(12) 



(13) 
(14) 

(15) 



In formulas (f8))- (|15l) <? MI/ is the metric tensor; _H" M " = + u^u v ; = d/(dx^); P, e, S, and /i c are the pressure, energy 
density, entropy density, and relativistic electron chemical potential, respectively. These quantities are related by the formula 



-e + n c n c + fii-rii + TS. 



(16) 



Finally, r\ is the shear viscosity coefficient and £i n , £2, £311, £4n are the bulk viscosity coefficients. Because of the Onsager 
symmetry principle, one has 



(17) 



Moreover, if the bulk viscosities are generated solely by the direct or modified URCA processes, one has an additional constraint 
( Gusakovll2007l ) 



£ln = ?2£3n. 



(18) 



In the absence of superfluidity the only nonzero coefficient is £2 - the ordinary bulk viscosity. 

To close the system describing superfluid hydrodynamics one should put two additional constraints on the four-vectors 
■u? and w 1 ? > , 

(n)> 



(19) 

(20) 



The first constraint is the standard normalization condition while the second one indicates that the com oving frame, in which 



we measure various t hermodynamic quantities (e.g.,ni,e, ...), is defined by the condition it M = (1,0,0,0) | Gusakov fc Andersson 



2006; lGusakovll2007h . Using Eqs. ([TJ, (lll|l . (|12[l . (I19p . and (|20p one then immediately finds that n; = — MjiJm Q =n > P> e ) an< l 



Making use of the hydrodynamics described above, one ca n derive the ent ropy generation equation, valid for superfluid 
matter. Following the derivation of the similar equation (33) in I Gusakov! ((2007J), one arrives at 



T 



T 



where the entropy density current is0 



= Su" - — T^. 
T 



(21) 



(22) 



Notice that, in iGusakovl (120071 1 there is an additional term in the expression for 5 M , so that 

The last term here appears naturally in the entropy generation equation. However, strictly speaking, it is small and should be neglected if 
one takes into acco unt only the largest diss ipative terms in the equations of superfluid hydrodynamics (this is the standard approximation; 
see iGusakov 200]] and §140 of lLandau &: Lifshitz 1987 for an explanation of what we mean by the 'largest terms'). It remains to note 
that the terms similar to the last term in the expression for also appe ar in the most ge neral form of the nonrelativistic superfluid 
dissipative hydrodynamics formulated by Clark (for details see the book by |Puttermanlll974|y 
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When writing (|21[) we neglected small dissipative terms, as it is discussed in Gusakov |2007l ). Introducing 



Qbuik = \V^\Y nk w? k) \ +V6<Mr . ( 23 ) 

Qshcar = TjH^ 1 H" 6 (^U-, ]S + U S]1 ~ | g-,6 U„ ;M 

= | H^ 1 H uS (u-y-s + u Sn - | fit„ ;M + u M; „ - | g M „ u%^\ , (24) 
Eq. (|21|) can be rewritten as 

T(S'u M ) ;fl = Qbulk + Qshcar- (25) 

To derive Eq. (f25j) we used Eqs. (fTT)) and (fT8)l . as well as the fact that for the tensor (fT2|) u„ =00. 



3 BASIC EQUATIONS 
3.1 An unperturbed star 

An eq uilibrium configuration of a nonrotating superfluid NS was analyzed in detail in section 3 of Gusakov fc Anderssonl 



I 2006). Here we present only the main results of this analysis, which will be used in what follows. 



The metric of a spherically symmetric, nonrotating NS in equilibrium has the form 

- ds 2 ee g [ °l&x a &x P = -edt 2 + e x dr 2 + r 2 (d6 2 + sm 2 8d<p 2 ), (26) 

where r, 8, and <p, are the spatial coordinates in the spherical frame with the origin at the stellar centre; t is the time 
coordinate; u(r) and A(r) are the metric coefficients for an unperturbed star. 
The four- velocity , generally defined as 

dx M 

u" = (27) 
ds 

in equilibrium equals 

-k/2 1 2 3 n /qq\ 

U = e , u — u = u — 0. (28) 

We assume that in the unperturbed star superfluid components are at rest with respect to the normal component. In that 
case the four-vectors satisfy 

W U = w fa = 0- (29) 
Using Eqs. Q, ([5]), (|28|l . and (|29|l . one has for the baryon four- velocity 

Ufa = e-" /2 , Ufa = Ufa = Ufa = 0. (30) 
In addition, the following conditions of hydrostatic equilibrium must hold for an unperturbed star, 
d-P 1 , „ , Au 

d7 = "2( p + e ^' ( 31 ) 

i {^' 2 ) = 0- (32) 
The last c ondition should be only use d in the stellar region where neutrons are superfluid (hereafter the SFL-region). One 



can show ( Gusakov fc Andersson 2006]), that if an unperturbed NS is additionally in beta-equilibrium, that is, the imbalance 



<5/z of chemical potentials vanishes, 

Sfi = /i n — fip — n c — 0, (33) 

then the SFL-region must also be in thermal equilibrium, with the redshifted internal stellar temperature T°° constant over 
this region, 

T°° ee Te v/2 = constant. (34) 

In what follows we assume that the conditions (|33[) and (|34[) are satisfied in the entire core of the unperturbed NS. In the 
latter case Eq. (|16[) for the equilibrium pressure can be rewritten as 

P = -e + ^ n n b + TS. (35) 

It should also be stressed that, as long as we neglected the temperature effects when calculating the equilibrium stellar model, 



2 This equality holds true only if one neglects the thermal conductivity, as we assume in Eq. 1121 . 
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the hydrostatic structure of the unperturbed superfluid NS is indistinguishable from that of the normal (nonsuperfiuid) star 
of the same mass. 



3.2 Small departures from equilibrium 

The metric of a perturbed star can be presented in the form 

- ds 2 = gapdx^dx 13 = (g^j + 6g a p)dx a dx 13 . (36) 

From here on the symbol 8 denotes Eulerian perturbations, so that 8g a p corresponds to small metric perturbations in the 
course of stellar oscillations. 

Since we study oscillations of a nonrotating nonmagnetized NS and neglect the effects of crystalline crust, all the pertur- 
bations in the system are of even parity n . In that case, in the appropriately chosen gauge 8g a p dx'^dx 13 can be written as (we 



follow the notations of Cutler et al. 199 1 



Sg a p dx a dx f3 = -e"r l H (r)Y l m e iu ' t dt 2 ~2iujr l+1 H 1 {r)Y t m e iu ' t dtdr 

-e A r l H 2 (r) e iu}t dr 2 - r l+2 K(r) Y, m e 1 "' (d(9 2 + sin 2 6d<p 2 ). (37) 

In Eq. (|37[) we assumed that all the perturbations depend on t as e 1 "'. In addition, we already expanded the perturbations 
into series in spherical harmonics Yj m , and consider a single harmonic with fixed / and m. The unknown functions Ho, Hi, 
H2, and K depend on r only, and should be determined from the linearized Einst ein equations, desc ribing NS oscillations (see 



Sec. [5j. Depending on I the gauge of the metric can be further specialized (e.g., Cutler et al. 1990t ). Namely, one can choose 
the gauge such that for I — (radial oscillations) Hi = K = 0; for / = 1 (dipole oscillations) K — 0; for / ^ 2 Ho = Hi. 

As follows from the definition (I27[) . in the perturbed star the four- velocity u M of the normal component equals, in the 
linear approximation 

u° = —! = =e-'/ 2 (l-lr l H Yre iut ), (38) 

V-goo V 2 j 

u j = v j e~ v/2 , (39) 
where 

v* S % (40) 

is the j-th component of the velocity of the normal liquid. Here and below j is the spatial index, j — r, 8, and if. Similarly, 
using Eqs. (|19[) and (|20|) one can show that for small deviations from equilibrium 

= 0, (41) 



while the spatial components w 3 ^ are small quantities, linear in perturbation (for a similar consideration see I Gusakov fc Andersson 



_ w ... . 

2006). In what follows, instead of the four- vectors [which are constrained by Eq. ©] it will be often more convenient to 



use the quantity A M , defined by @. For small perturbations 

X° = 0, (42) 

while X 3 is non-zero but small (linear in perturbations). 

Using Eqs. (|38p . (|39|l and (|42[). as well as the definition ©, it is easy to write out an expression for the baryon four- velocity 
Ufo in the perturbed star, 

C/ ( °b) = - 7 L==e- v/2 (l-\r l HoYre 1 " t ), (43) 



V-ffoo \ 2 

'(b) - "(b) ' 



UL = vLe-" /2 , (44) 



where the last equality is the definition of the j-th component of the baryon velocity v J ^ (linear in perturbation). Notice 
that, as follows from Eqs. (|43[) and (|44|l . in the linear approximation the normalization condition for the baryon four- velocity 
is the same 

U Wl >Ufa=-l, (45) 
as for m m 0. In what follows instead of the velocities v J and vLs , it will be more convenient to use the corresponding Lagrangian 

3 A more detailed argument can be found in lThorne fc Campolattarol Jl967f) ; see also iRegge fc Wheeled Jl957t> . 

4 However, beyond the linear approximation, Eqs. Ht . JSJ, 1191 . and 1201 1 yield ^(b)^^(b) = — I + Y n iY n ^ ^w^/n^ . The nor- 
malization condition 14511 is generally not fulfilled because the reference frame in which = (1,0,0,0) is not comoving [that is, 
3(b)^0>) H ^ ~ n b in that reference frame]. As it was already indicated in Sec. [2] all thermodynamic variables are measured in the 
reference frame, in which n M = (1, 0, 0, 0). 
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displacements. They are defined by the equalities 

^=^1=1^, (46) 

w (*) s T =i <)- (47) 
Introducing also the analogue of the Lagrangian displacement & s m for the vector X 3 , one can write 

X^e-^^=i W e-^U- (48) 
In terms of the Lagrangian displacements the equality (JSJ can be presented as 

Because of the spherical symmetry of the unperturbed star it is sufficient to consider Lagrangian displacements £ J , 
and ^ sfl j, of the form [see also a note after Eq. (|90p below] 

e = [e,Z e ,e} = [e~ X/2 r l - 1 W(r)Y l °, -r'- 2 V(r) OeY?, O^e 1 ^, (50) 

e w = ffb)] = [e~ V2 ^ W,(r) if, -r 1 " 2 V h (r) fcY, , o] e iwt , (51) 

^ sfl) = [C( sfl) , = [ e ~ V2 ^ y '°' v « W °] ^ ( 52 ) 

where W, V, Wb, Vt>, H4a> and T4a are some functions of r to be derived from oscillation equations. In Eqs. (150[) (|52|l 



y ; ° = ^/ (2/ + 1)/ (47r) P; (cos 6) , where P; is the Legendre polynomial. Here and below we consider only spherical harmonics 
with m = 0. We can do this without any loss of generality, because, due to the spherical symmetry of the unperturbed star, 
oscillation ei genfrequencies as well as eigen f unctions Ho, Hi,. . ., W s a, and V s a, introduced in this section, cannot depend on 
m (see, e.g.. lThorne fc Camporattardll967l ). 

It follows from Eqs. (J49J) and M-EIl that 



W h = W + W s n, (53) 
V b = V + V s& . (54) 



4 DAMPING OF OSCILLATIONS DUE TO THE BULK AND SHEAR VISCOSITIES: GENERAL 
FORMULAS 

In the present paper among the possible mechanisms of dissipation of oscillation energy we take into account damping due to 
the bulk and shear viscosities as well as due to radiation of gravitational waves. Dissipation makes the oscillation frequency 
cj complex, so that it can be presented in the form, 

Lj = a + -, (55) 

T 

where a is the real part of the frequency, and t is the characteristic damping time. Assuming that damping is weak, in the 
linear approximation one can present the following standard expression for r, 

1 1 d_B mcc h , . 

7 = ~9P *i~> ( 56 > 

t /_t/ mcch at 

where _E mcc h is the mechanical energy of oscillations; dP mec h/d£ is the dissipation rate of the mechanical energy, which can 
be presented as 

dE ™ Ch = -2Ubulk - 2IWar - 2U grav , (57) 

where 2Bbulk, 2Ushcar, and 2U gra v are the energy, dissipated per unit time due to the bulk viscosity, shear viscosity, and 
gravitational radiation, respectively. Introducing partial damping times Tbuik, T s hear, and r gr av according to 

1 = gPbulk 

1 2Ughcar 



Tshcar 2£/ me ch 
1 2H sr 



7~grav 2_£/ m cch 

one can rewrite the expression for r as 



(58) 
(59) 
(60) 



I = J- + _^ + l (61) 

» 7~grav ^shcar ^"bulk 
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Thus, to calculate r we need to know the mechanical energy -E mccn of NS oscillations, as well as the quantities 2Ubuik, 2H s hcar, 
and 23grav 

4.1 Mechanical energy 

The general relativistic expres sion for the mechanical energy o f oscillating normal (nonsuperfluid) NS was obtained by 
Thorne fc Campolattaro (1967) (see also Meltzer k, Thorne 1966T I. Their result can be easily generalized to the case of super- 
fluid matter. Mechanical energy E mccn is related to the averaged over the oscillation period 2-k/g kinetic energy -Ekin by the 
standard formula, 



-En 



= 2 -E kill 



Thus, to determine -E moc h one needs to know -Ekin- One can write (e.g., iThorne fc Campolattardll967r ) 



-Ekin = / Ekin e 
' star 



■12 



dV, 



(62) 



(63) 



where dV = r 2 e x ^ 2 sin#d#d<y3dr is the proper volume element; ekin is the kinetic energy density measured in the locally flat 
coordinate system x^ [with the metric — ds 2 = g^dx^dx" , where (? M „ = diag(— 1, 1, 1, 1)], which is at rest with respect to 
the unperturbed star. If the NS matter is normal, ekin is given by 



ekin = ^(P + e) \(uy + (u°f + (un 2 



Here 



dx 3 
dx^ ' 



-v/2 A/2 r 9 

e e ' v , r v 



r sin# v* 



(64) 



(65) 



is the physical velocity of the fluid in the locally flat coordinate system x^ . For superfluid matter Eq. (|64[1 should be modified, 
because in this case not only m otion of the normal liquid component contribute to eki n but also that of the superfluid 
component. Using formula (56) of Kantor fc Gusakov ( 20091 ) one obtains [f] 



ekin = 



- |(P + £)[m j ] 2 + Y ik [tMW 3 ik) Uj +fj,k w^Uj +i5( i) *(fc)j] } ■ 



where w 3 ^ = [dx 3 /dx^w^. Taking into account Eqs. ((3}-([5]) and (|35[1 . Eq. (|66|) can be rewritten as 



+ y 

where we neglected 'temperature' term TS in the expression (|35[) . In Eq. ((67 

nb^pp 



n 3 
u (b) 



x 3 



^n(^nn^pp ^np) 
9x 3 „ -v/2 f A/2 r ■ a in I 

— C/f b) =e |e v (b) , rv w , rsmSv^ 
dx 3 



dx» 



X p 



[< 



A/2 



X r , rX e , rsmeX" 



(66) 
(67) 

(68) 

(69) 
(70) 



Now, using Eqs. @7), {48]), |5T|I, and ([52j) let us express (J69J) and ((70]) through the functions W b (r), V h (r), W Bfl (r), and V B a(r), 
and then substitute Eq. (1671) for ekin into (|63|) . After integrating Eq. (|63p over sinOdOdip (in the same way as it was done in 
Thorne fc Campolattarolll967h and making use of Eq. (|62p . one arrives at the following expression for -E mec h, 

Smcch(i) = E? mcch ( b )(t) + E moch ( sfl ) (t) , (71) 

where we tentatively presented £? me ch as a sum of two terms related to the baryon motion as a whole -E mcch ( b ) and an 
additional term -E mech ( sf i) appearing because of the superfluid motion, 

Ft 



E 



mcch (b) 



(*) 



2 -2t/r 

■ o e 



(P + e) e 



(\-v)/2 21 



r 21 [W^ + 1(1 + 1) V b 2 } dr, 



E, 



mcch (sfl) (i) — g ° 6 



- 2t/T [ R (P + e)e 
Jo 



(A-^)/2 21 



r 21 y [W 2 a + 1(1 + 1) Kl] dr. 



(72) 



(73) 



Strictly speaking, the functions W b (r), V b (r), W s fi(r), and V$fi(r) in these formulas are complex, that is, instead of, for example, 
W b (r) 2 one should write | WbMI 2 - Notice, however, that all these functions [as well as Hq(t), Hi(r), H2(r), and K(r)] are 
defined up to the same arbitrary complex multiplicative constant. Since a>l/r (dissipation is weak), one can always choose 
the constant in such a way, that the real parts of all these functions would be much greater than their imaginary parts (e.g., 



5 This exp ression is analogous to th e corresponding formula for the kinetic energy density of a nonrelativistic superfluid mixture, 
obtained bv lAndreev &z Bashkin (1975|), see their equation (7). 
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He[H2{r)] 3> lm[H2 (r)]), so that one could neglect their 'complexity'. From here on, unless otherwise stated, by the functions 
Wb(r), Vb{r), W s fl(r), V B n(r), Ho(r), Hi(r), H2(r), and K(r) we mean their real parts. 

In the absence of superfluidity W b q = V s a = 0, Wb = W, and Vb = V. In that case Eq. (|72l) g ives a mechanical energy 
of a n onsuperfluid star that coincides, up to notations, with the corresponding expression (29) of Thorne fc Campolattarol 
I l967f >. 



4.2 Dissipation rates 



The damping time r grav due to radiation of gravitational waves can be obtained from the equations, describing linear oscilla- 
tions of NSs (see Sec. [5] below). The goal of the present section is to determine the dissipation rate of oscillation energy due 
to the bulk 2H bu ik and shear 2U s hcar viscosities and, as a consequence, the damping times Tbuik and r s h car . 

For that, we turn to the entropy generation equation (|25[1 . Using it, one can find rate of change of the (averaged over the 
oscillation period) thermal en e rgy o f a star dE t h/dt due to bulk and shear viscosities. Following the derivation of Eq. (34) in 
Gusakov. Yakovlev fc Gnedinl (120051 ). one obtains 

dE th 



dt 



(74) 



where Q bulk and Q shear are the values of Qbuik and Q s hear, averaged over the oscillation period 2ir /a [see Eqs. (|23p and (|24|) ]. 
Obviously, the increase in the thermal energy E t b is accompanied by the decrease of the oscillation energy ,E mcc h, that is 



2B bu i k = / Q buIk e"dy, 



an, 



shear — / Vshear 
star 



Qshear ^V. 



(75) 
(76) 



Using these equations, as well as the formulas (|23l) . (|24[1 . I|58[). (J59J) , and the definitions of Sec. 13.21 one gets, after rather 
lengthy calculations, 



Tbulk 
1 

7"shcar 



where 

AM 

ai(r) 
a 2 (r) 



4S m ech(0) 
2£mech(0) Jq 

X || (aif + 21(1 + 1) (a 2 f + 1(1 + 1) 



dr, 



2(1-1) A/2 

r\ r v ' e ' 



V > dr, 



f i 1 r_r 1 -A/2 

K + -H 2 e 1 

2 r 



dW 1,, 
dr r 



1 

r 


-A/2 

e ' 


r 2 




3 




r 


"dU 


2 


dr 



d(n b W / s fl) 1,, . ' 

'- + - (I + 1) n b Wsh 

dr r 



1(1 + 1) 



A/2 



^ + (1-2)^' 
dr r 



2 ' 
n h U s fl 



y 

r 



w 

r 



+ K - H 2 - 1(1 + 1) ^5 

-A/2 



(77) 

(78) 

(79) 
(80) 
(81) 
(82) 



As for the mechanical energy (|71[) . to obtain from these formulas r bu i k and r s hear for a nonsuperfluid star, on e has to put 



Wad = VsB = 0. In that case our Eqs. (|77[) and (|78|l should coincide with the corresponding formulas (5) and (6) of lCutler et al 



(1990). Unfortunately, direct comparison of these formulas reveals, that our r bu ik and Tshear appear to be 2 times larger 
Using, as tests examples, damping of: (i) NS radial oscillations, (ii) p-modes i n the NS envelopes, and (Hi) sound waves in 
the nonsuperfluid matter of NSs w e checked, that our results reproduce those of Gusakov et al. ( 20051 ): Chueunov fc Yakovlev 
( 2005h : lKantor fc Gusakol |2009h . obtained in a quite a different way. 



5 OSCILLATION EQUATIONS 



In order to calculate r bu ik and Tshear one has to determine the oscillation eigenfrequencies a and eigenfunctions Ho, Hi, H2, 
K, Wb, Vb, W s fi, and V s a. To do that one needs to formulate oscillation equations. Since the dissipation is weak, when deriving 
the oscillation equations one can neglect the dissipative terms in the superfluid hydrodynamics of Sec. [2] and put = and 
x n = 0. 

As it was shown in lGusakov fc Kantorl (|2011i ). equations, describing small linear oscillations of a NS include: 
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(i) Continuity equations for baryons |[5J and electrons ([9)l , that can be written in terms of the baryon and electron number 
density perturbations, <5rib and Sn e , as 



(83) 
(84) 

(85) 
(86) 



Snb = Z^^[ dj{nh)U ^ +nhU ^^}' 

8n c — 5n e (norm) + (sfl) 5 

where j is the spatial index and we defined 
5n c (norm) = ^ e ^ /2 \dj (n c ) U J (b) + n c Ufa ;M 

5n e ( sfl ) = - \dj (lie) X 3 + n e 

(ii) Einstein equations, which can schematically be presented as 

6(1^" - 1/2 g^ R) = 8ttG 5T^ , (87) 

where the perturbation ST 11 " of the energy-momentum tensor (JTTJ can be expressed in terms of the perturbations of baryon 
four- velocity SUfa, metric <5g M „, pressure SP and energy density Se as 

ST" U = (5P + Se) UfaUfa + (P + e) [ufa SUfa + Ufa SUfa] + SP + P Sg^. (88) 

In Eq. (187[l R^" and R are the Ricci tensor and scalar curvature, respectively; G is the gravitation constant. 

(Hi) 'Superfluid' equation, that can be derived from Eqs. (|10p and (| 13 j) of Sec. [2] (here we present only the spatial 



components j of this equation) |f| 



ioj (fi n Y nk W(j,)j — n b iW(n)j) = n c 9j(e^ /2 <5/i). 



(89) 



Expressing the vectors w 3 ^ through X 3 in this equation [see Eqs. and Q], and introducing the redshifted imbalance of 
chemical potentials 5fi° 



Sfi, one can rewrite Eq. (|89p as 

i n c 



■9j(Sn° 



(90) 

/x n n b w y 

where y is defined by Eq. (|68[) . Notice, that this equation dictates the most general form of the superfluid Lagrangian 
displacement & s m, that was already obtained in Eq. (|52|) from the symmetry arguments. 

Eqs. ()83l)- (|90[) should be supplemented with the expressions for the perturbations 8P, 5fi, and 8e. To derive them, let us 
notice that any the rmodynamic q uantity (e.g., P) in the superfluid matter can be presented as a function of rib, n e , T, and 
Wfji ./W^-) (see, e.g., Giisakov|l2007 ) . In strongly degenerate matter the dependence of P, Sfi, and e on T can be neglected (see, 
e.g.. lReiseneggerlll995l : iGusakov et al.ll2005l ). while the scalars w^^wfa are quadratically small in a slightly perturbed star 



[see Sec. 13. 2j . Thus, P — P(n^, n c ), <5/i = <5/i(nb, n c ), and e = e(n^, n c ). Expanding these functions into Taylor series near 
the equilibrium, one obtains 



8P = 

8/i — 

Se = 



rib 



dP 



dn h 

dSu, 

n <= a — 
an c 

/Xn Srib- 



Sn h 



rib 



Sn 



c (norm) 



+ S ■ 



Sn, 



(sfl) 



Sn 



e (norm) 



+ 



Sn, 



o (sfi) 



He 



where we made use of Eq. 



and introduced dimensionless coupling parameter s and the quantities s and z, 



(91) 

(92) 
(93) 

(94) 
(95) 
(96) 

Notice that the variable s is equal to s here. The reason for discriminating between s and s is purely technical: To solve 
oscillation equations (see Sees. [6] and it turns out to be convenient to develop a perturbation theory in (small) parameter 
s, at the same time treating the terms depending on s in a non-perturbative way (see Sec. 16.21 and, in particular, footnote^ 



n e 


(dP/dn c ) 


n b 


(8P/dn h ) ' 


n c 


(dP/0n e ) 


n h 


(8P/dn h ) ' 


n h 


(d5fi/dn h ) 


n c 


(dSfi/dn c ) 



6 It is worth to make a number of comments on Eq. l[89ll : (i) In lGusakov fc Kantoj \20lt \ this equation was derived under the assumption 
that the only s uperfluid species are neutrons (that is Y p i = 0). A generalization of this equation to the case o f poss i ble proton superfluidity 
is pre sented in lChugunov fc Gusakovl i2011f ) [see their Eq. (3)]; (ii) In both papers. [Gu sakov fc Kanto"n i2011f ); [Chugunov fc Gusakovl 
(2011), this equation is written with the same mistake. In particular, in IChugunov fc Gusakovl l|201ll l one should write n c dj (e"/ 2 Sfi) 
instead of rice"/ 2 dj(8fi) in the right-hand side of Eq. (3). 
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there). When deriving Eq. (|93l) we took into account that [<9e(nb, n e )/dn e ]Sn c — —Sjj,Sn c is a quadratically small quantity, 
because Sfj, = in equilibrium [j. 

The vector superfluid equation (|90p can be substantially simplified, and reduced to a scalar one. For that let us notice 
that, without any loss of generality, the scalar <5^i°° can be presented as 

5n°° = 6 f ii(r)Y l °(e). (97) 

Employing now Eqs. (|90|) and (|92[) . one arrives at 



<w + 



A' 2 



1(1+1) -„/a U 

r 2 h«8 



/iQ3 



(98) 



Here ft = e"/ 2 n 2 /(/i n rib y), 23 = d8[i(nh,n e )/dn a , and prime (/) means derivative with respect to the radial coordinate r. 
Furthermore, <5/x n orm i (r) in Eq. ((98} is defined by 

<5/-Crm = <W>rmi( r ) (99) 

where 



i//2 



5n b , 

2 h 

n b 



6n, 



c (norm) 



(100) 



is a part of 5fi°°, which depends on &g^ v and and is independent of the superfluid degrees of freedom X 3 [see Eqs. (|83|) 
and (|85[) ]. The function (fyinormiM can be easily rewritten in terms of Ho(r), Hi(r), H2(r), K(r), Wb(r), and Vb(r) with the 
help of Eqs. (|33), (gll, HI), (gll), HU, J83]), and (|85}. One obtains 



A// 



„/ 2 dS[i(n b ,x c ) i 

normi = e U b T fix, 



(101) 



where x c = n c /nb and /3i(r) is given by Eq. ()79[) with Wb and Vb instead of, respectively, W and V. 

Finally, let us mention one important property, that follows from the oscillation equations and quasineutrality condition 
|(3}. If neutrons in a nonrotating nonmagnetized star are normal (i.e. Ynn = Y np — 0), while protons are superfluid (Y pp 7^ 0), 
then oscillation eigenfrequencies and eigenfunctions for such star will be indistinguishable from that for a normal star of the 
same mass (where both protons and neutrons are nonsuperfluid). 



6 OUR APPROACH 



6.1 Decoupling of superfluid and normal modes 



In principle, Eqs. (|83[) - ()101[) allow one to study the nonradial oscillations of superfluid NSs and thus to determine the 
spectrum of eigenfrequencies uj, eigenfunctions Ho, Hi,. . ., W s n_, and I4g, and hence the damping times T gr av, Tbuik, and r s h ear . 
However, this task can be significantly simplified, i f one notes that the dim ensionless coupling parameter s ()94|) is small 
for realistic equations of state of superden se matter (jGusakov &i Kantorll201lT ). For example, for the equation of state APR 
l Akmal. Pandharipande fc Ravenhall 19981 ) employed below s ~ 0.01 -f - 0.05. This means that one can look for the solution to 
the system of Eqs. ()83[) (IXOTT ) in the form of a series in s. Since s is small, the approximation s = is already quite accurate. 
Indeed, as it was shown in lGusakov fc Kantorl (|2011f ) with the example of radial oscillations, the eigenfrequencies calculated in 
this approximation differ from the exact ones, on average, by ~ 1.5-^2%. Thus, in what follows all calculations are performed 

assuming s — 0. 

How this approach simplifies the problem? As it was first demonstrated in iGusakov fc Kantorl (|20111 ). in the s — 
approximation superfluid degrees of freedom (vectors X 3 ) completely decouple from the 'normal' degrees of freedom [metric 
perturbations Sg^v and baryon four-velocities SUj^.]. That is, one has two distinct classes of oscillations: 'superfluid' and 
'normal' modes, which are described by independent equations. For superfluid-type oscillations the metric and baryon velocity 
are not perturbed [dg^ = and SU^ = 0] , hence these modes do not emit gravitational waves; moreover, they are entirely 
localized in the SFL-region. At the same time, the frequencies of normal modes are indistinguishable from those of a normal 
(nonsuperfluid) star of the same mass[f|. Below we discuss in more detail decoupling of superfluid and normal oscillation modes 
and how this property can be used to calculate the characteristic damping times. 



7 The equality 9e(rib, n c )/dn c = —5fl follows from the second law of thermodynamics H15I) . which can be rewritten in our case as 
de = u n drib — "5 A 1 dn e . 

8 Here and below by 'normal modes' we mean oscillation modes (of approximate solution), that also exist in the normal (nonsuperfluid) 
star. 
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6.2 A strategy to calculate the damping times 

So, let us formally assume that s = (while s is given by Eq. (|95|l and is non-zero). Then, as follows from Eq. (|91ft . SP equals 

(102) 



dP 

5P = n b ^- 



8n h _ Sn e ( norm ) 

Mb Tie 



and is independent of the superfluid degrees of freedom X th [see Eqs. (|83ft and (|85|> ] . Other terms in the expression (|88|l for 
(jT'"' also do not depend on X M [in particular, 8e = fi n Snb does not depend on X^ due to Eqs. (|83ft and (193ft ] . Thus, we 
come to conclusion that the linearized Einstein equations (|87ft depend only on perturbations of the metric g M „ and the baryon 
four- velocity U^ h ) and are independent of X 1 *. Moreover, it is easy to see, that in the case s = these equations (and the 
corresponding boundary conditions) have exactly the same form as in the absence of superfluidity 0. Correspondingly, two 
alternatives are possible when solving the system of Eqs. (|83ft - (|101ft in the approximation s = 0: 

(1) A star oscillates at a frequency which is not an eigenfrequency of the Einstein equations 1)871) . In that case, to satisfy 
Eq. (|87ft . one has to demand 

Ho = Hi = H 2 = K = W b = V b = 0. (103) 
From Eq. (|10ip it follows then, that i/i norm i = and the superfluid equation (|98f) decouples from t he Einstein equations 



As a r esult we arrive at the 'source- free' equation (with the right-hand side vanished), first derived in lChugunov fc Gusakov 

B0 



'V'" - ( \ - \ + \ ) 'V'.; - ■ 1 



r 2 h<B 



6 fit = 0. (104) 



This equation describes superfluid modes and should be solved in the stellar region where neutrons are sup e rfluid (SFL-region 



It should be supplemented with a number of boundary conditions, discussed in IChugunov fc Gusakovl (|201lf ) [similar, but 
more general boundary conditions for Eq. (|98p are presented in the Appendix]. Having solved Eq. (|104p for Sfii, it is easy to 
determine the functions W a a and Vfl using Eqs. (|48p . (I52p . and (I90p . Using W sn and V sn , one can find the functions W and 
V from Eqs. ([53]), ((54j) , and (IT03|l . 

W=-W s n, V = -V sfl . (105) 

This information is sufficient to calculate rbuik and T s h car from Eqs. (177ft and (|78p [as follows from Eq. (|103p . r grav = oo for 
superfluid modes in the s = approximation]. 

(2) A star oscillates at a frequency which is an eigenfrequency of Einstein equations (|87p . In that case, the eigenfrequency 
and eigenfunctions Ho, Hi, H2, K, Wb, and Vb are indistinguishable from the corresponding eigenfrequency and eigenfunctions 
for an oscillating nonsuperfluid NS [we recall, that for the nonsuperfluid star Wb = W , Vb = V , because W s fi = V s a = 0, see 
Eqs. (|53p and (|54ft]. There is, however, one very important difference: for a superfluid star the functions W a a and V a R do not 
vanish in the SFL-region and are comparable there to Wb and Vb. As follows from Eqs. (|77p and (|78ft . the damping times 
Tbuik and T s hoar depend on these functions [as well as on W = Wb — W a a and V = Vb — V s r], that is why the determination of 
Wafi and Vfl is a necessary task. 

To determine these functions we make use of Eq. (|98ft . Since the oscillation frequency uj — er+i/rgrav and the eigenfunctions 
Ho, Hi, H2, K, Wb, and Vb are already known, we can, using Eq. (|101ft . calculate i/i, lorm i and determine a 'source' in the 
right-hand side of Eq. (|98ft . This source plays a role of an external driving force, that makes the superfluid equation (|98ft 
'oscillate' at the frequency uj, which is not an eigenfrequency for this equation Lj. As a result, the function 8fii(r) will be 
nonzero. To determine it one has to specify the boundary conditions for Eq. (|98ft ; they are formulated in Appendix. Having 
solved Eq. (|98ft numerically and having defined Sfii(r), one can calculate the functions W a a and V s a, using Eqs. (|48ft . (|52l) . 
(|90ll. and (|97|). 

Summarizing, in the approximation s = the eigenfrequencies and eigenfunctions Ho, Hi, H2, K, Wb, and Vb, (and 
hence r grav ) for the normal modes appear to be the same as for a nonsuperfluid star. At the same time the eigenfunctions W a a 
and V s fl are non-zero in the SFL-region and should be determined from Eq. (f98)l . As a result, the damping times Tbuik and 



9 This is the main advantage of treating s in a non-perturbative way. Notice, however, that this trick leads to somewhat 'excessive' 
accuracy of the approximate solution to oscillation equations: the retained terms depending on s may lead to smaller correction to 
the solution than the ^-d ependent terms which were ignored. Bearing this in mind and with the aim to simplify consideration, in 
iGusakov fc Kantorl {2011) it was assumed that both parameters s and s vanish in the s=0 approximation. Such an approach is also 
possible. In that case, strictly speaking, the resulting Einstein equations would differ slightly from the equations describing oscillations 
of a nonsuperfluid NSs. In particular, instead of the standard adiabatic index of the 'frozen' npe-matter 7f r = (nb/ P)[dP{nb, x c )/dnb] 
the new index would appear, 7 = (nb / P)\dP(rib,n e )/dnb]- Howeve r, this difference is not essential, because s and s are small. 

10 The corresponding equation (5) of Chugunov & Gusakov (2011), contains a mistake, that was corrected in the second version of the 
manuscript in arXiv (see arXiv:1107.4242v2). 

11 In the present paper, in all numerical calculations we used <r instead of u in Eq. l|98jl . because a 3> 1/rgrav- Also, when calculating 
8fi nornl [ we only employed the real parts of eigenfunctions Ho, Hi, Hi, K, Wb, and Vb [see a note after Eq. 1731 1 . 
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Figure 1. (color online) Left panel: Nucleon critical temperatures T c k (k = n, p) versus density p for model 1. Right panel: Redshifted 
critical temperatures T£g versus radial coordinate r (in units of R) for model 1. 



Tahear, defined by Eqs. (|77[) and (|78[) , will differ from the corresponding times, calculated using the ordinary (nonsuperfluid) 
hydrodynamics (even if one takes into account the effects of superfluidity on the kinetic coefficients). 



7 RESULTS 

Let us apply the approach, suggested in the previous section, to determine the frequency spectrum and damping times for an 
oscillating superfluid NS. But first let us discuss its equilibrium model. 



7.1 Microphysics input and equilibrium model 

As m entioned in Sec. El we consider the simplest npe-com position of NS core. We adopt APR equ ation of state J Akmal et aL 
199ct ) parametrized by Heiselberg fc Hiorth-Jense nf j 19991 ) in the core and the equation of state bv lNegele fc Vautherinl (|l973l ) 
in the crust. 

All numerical results presented here are obtained for a NS with the mass M = IAMq. The circumferential radius for such 



from the centre. 

When modeling the effects of superfluidity we assume the triplet pairing of neutrons and singlet pairing of protons in the 
NS core. The neutron superfluidity in the stellar crust is neglected; it should not affect strongly the global oscillations of NSs. 

We consider two models of nucleon superfluidity: model '1' (simplified) and model '2' (more realistic). In the model 1 the 
redshifted proton critical temperature is constant over the core, = T cp e" /2 = 5 x 10 9 K; the redshifted neutron critical 
temperature = T cn e"'' 2 increases with the density p and reaches the maximu m value T™ max = 6 x 10 8 K at the stellar 
centre (r=0). This model corresponds to the model 3 of iKantor fc Gusakovl l|201ll ). 

In the model 2 both critical tem peratures T cn and T cp are density dependent. This model does not contradict the results 



of microscopic calculations (see, e.g., Lombardo Sz Schulzel 200ll : Yakovlev et al. 19991 ) an d is similar to the nu cleon pairing 



models used to explain observations of the cooling NS in Cassiopea A supernova remnant (jShternin et al.ll201ll ). 

The models 1 and 2 are shown in Figs. [1] and [2] respectively Q The function T C i(p) in both figures is shown in the 
left panels, while the right panels demonstrate the dependence T£°(r) [i — n and p]. With the decrease of the redshifted 
temperature T" 30 the size of the SFL-region [given by the condition T < T cn (r), or, equivalently, T°° < T^(r)] increases or 
remains unchanged. For instance, the SFL-region corresponding to T°° = 4 x 10 8 K, is shaded in Figs. \T\ and [2] One can 
see that for the model 2 there can be three-layer configurations of a star with no neutron superfluidity in the centre and in 



Fig. [2] is a slightly modified version of Fig. 1 from lChugunov fc Gusakovl ll201ll) . 
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Figure 2. (color online) The same as in Fig. [T]but for model 2. 



the outer region but with superfiuid intermediate region. On the contrary, in the model 1 only two-layer configurations are 
possible. 

The entrainment ma trix is calculated for the superfluidity models 1 and 2 in a way similar to how it was done in 
Kantor fc Gusakovl (|201lh . 

When analyzing viscous dissipation in oscillating NSs we allow for the damp ing due to shear and bulk vi scosities. For the 
shear viscosity coefficient r\ we take the electron shear viscosity rj c , calculated in lShternin fc Yakovlevl (|2008l ) . We neglect the 
nucleon shear viscosity because: (i) it is poorly known even for nonsuperfl uid matter and (ii) it appears to be less than the 

electron shear viscosity in the core at T « T cp [jShternin fc YakovlevlExil). 

The bulk viscosity coefficients are calculated as described bv lGusakovl (|2007l ); Gusakov fc Kantor ( 20081 ) ; Kantor &: Gusakov 

(2011). Since the direct URCA process is closed for our stellar model with M = IAMq, the main contributor to the bulk 
viscosity is the modified URCA process. 



7.2 Oscillations of a nonsuperfluid star 



As follows from Sec. 16.21 before considering oscillations of a superfiuid NS one should study those of a normal (nonsuperfluid) 
star of the same mass. To this aim, we have determined the eigenfrequencies and eigenfunctions of the radial and nonradial 
oscillation modes for a nonsuperfluid NS of mass M = 1.4M a and equation of state APR (see Sec. I7.1|l . We have solved 
the equations describing radial and nonradial perturbations of a nonrotating star in general relativity. These equations are 
derived by expanding the perturbed Einstein's equations in tensorial spherical harmonics in an appropriate gauge, and are 
integrated in the frequency domain. 

Stellar modes are defined as solutions of the perturbed equations which are regular at the centre and with vanishing 
Lagrangian pressure perturbation at the surface, and (if I > 1) which behave as a pure outgoing wave at infinity; as discussed 
above, such solutions have complex frequencies uj — a + i/r. If I ^ 1, instead, the frequency is real and the mode is not 
associated to gravitational emission. 

The oscillation modes are classified according to the source of the restoring force which prevails in bringing the perturbed 
element of fluid back to the equilibrium position; for instance, we have a gr-mode if the restoring force is mainly provided by 
buoyancy, a p-mode if it is due to a gradient of pressur e, and so on. 



The radial modes are calculated as described in Gusakov et al. d2005l) . To calculate the nonradial modes we follow 



the formulation of Lindblom fc Detweiler (1983) and Detweiler fc Lindblom ( 19851 ). In their formulation, the equations for 
nonradial perturbations can be expressed, inside the star, as a system of first-order differential equations in the variables 
Ho, Hi, H2, K, Wh, and Vt> defined in Sec. [5] Outside the star, they reduce to a simple, second-order differential equation 
(t he Zerilli equation) . By numerical integration of these equations (the procedure we have followed is described in detail, e.g., 
in lBurgio et alj|201ll ) we find, for each value of the multipolarity I, the (complex) eigenfrequencies w and the corresponding 
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Table 1. Frequency a (in units of 10 4 s 1 and in units of <r = cj R PS 2.46 X 10 4 s x ) and the damping time T grav (in seconds) for various 
oscillation modes of a nonsupcrfluid NS. The first column shows the multipolarity £ of modes and their names. 
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Figure 3. The function <S/^ norm j (in units of 1 7 kelvins) versus r for fundamental radial F-mode as well as for pi- and /-modes with 
multipolarities I = 1, 2, and 3 (see the footnote . The energy of each oscillation mode is 10 43 erg. Shaded region corresponds to crust, 
where Sfi norm i is not defined and was not plotted. 

eigenfunctions Ho(r), Hi(r), H2(r), K(r), Wb(r), and Vb(r). The results of our computations are summarized in Table [1] and 
illustrated in Figs. [3] and [4] 

Table 1 presents the real parts of the eigenfrequencies Re(cj) = a (measured in units of 10 4 s _1 and in units of a = c/R « 
2.46 x 10 4 s _1 ) and the characteristic gravitational damping times r grav (in seconds) for the modes with I = (fundamental 
F-mode and first two overtones 1 H and 2 H), I — 1 (dipole pi-mode), I — 2 (quadrupole /- and pi-modes), and I = 3 (octupole 
/- and pi-modes) One can see, that a 3> 1/T gra v in all these cases. That is, damping due to emission of gravitational waves 
occurs on a time scale much longer than the oscillation period. 

Using the definition (|99[) and Eq. (|101[) we have determined, in terms of the eigenfunctions Ho(r), Hi(r), . . . , Vh(r), the 
function 5/x n0 rm i (r) and, consequently, the quantity <5/i r ^ rm (r, 6) = Sfj, norrn i(r) Yp(6) for each mode. As follows from Eqs. (|92[) 
and (|100|) . for a nonsuperfluid star 8fi^ rm (r, 9) is simply a redshifted imbalance of chemical potentials, S/j,°° = 5/^ rm . The 
function <5/i norm ;(r), entering Eq. (|98|) . is shown in Fig. [3] for the oscillation modes from Table 1. It is normalized such that the 
mechanical energy of oscillations is 10 43 erg. The shaded region corresponds to the crust of the star, where 8fj, norm i(r) is not 
defined (protons are bound in nuclei there). As seen in the figure, \Sfi nolln i(r)\ for /-modes is about one order of magnitude 



The /-mode is absent in case of I = 1. 
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Figure 4. (color online) Damping times Tb+ S = (tJ~Vl + r s ~hoar) _ 1 versus T 00 for various oscillation modes. The effects of superfluidity 
are partially taken into account, as described in the text. Thick and thin solid lines correspond to radial (I = 0) F- and 1 //-modes, 
respectively; dot-dashed line - to dipole (I = 1) pi-mode; thick and thin dashes - to quadrupole (I = 2) /- and pi-modes, respectively; 
thick and thin dots - to octupole (I = 3) /- and pi-modes, respectively. 



smaller than for p-modes. This is not surprising, since matter is only weakly compressed during /-mode oscillations, so that 
a deviation from beta-equilibrium (when 8fi°° = 0) is small. The functions 5/i llorm ; (r) are employed to calculate the damping 
times of a superfiuid NS in Sec. 17.41 

Fig. [4] shows the viscous damping time t^ +b = (t^^ + i^car) -1 as a function of T°° for a set of oscillation modes. 
The solid lines corresponds to radial (/ = 0) modes F and 1H; dot-dashed line to dipole (I — 1) mode pi; dashed lines to 
quadrupole (I — 2) modes / and pi; dotted lines to octupole (I — 3) modes / and pi. To calculate Tb+ S we used the formulas 
for Tbuik and r s h C ar, applicable for the ordinary hydrodynamics of a nonsuperfluid liquid0- However, we allow for the effects of 
superfluidity when calculating the kinetic coefficients rj and £2 (the other bulk viscous coefficients do not appear in the normal 
fluid hydrodynamics). To calculate r\ and £2 we adopt the nucleon superfluidity model 2 (see Sec. 17. Such an approximate 
approach to accounting for the effects of superfluidity is commonly used in the literature, but it is not fully consistent. The 
results of a more consistent approach (see Sec. [6| are discussed below in Sec. 17.41 

As follows from Fig. [4] the dependence of Tb+ a on T°° is a power-law at T°° < 6 x 10 8 K. At such T°° the proton super- 
fluidi ty is 'strong' (T°° <c T^). In that case the bulk viscosity is expon entially suppressed ijHaensel. Levenfish fc Yakovlev 
2001), while the shear viscosity r\ oc 1/(T°°) 2 jshternin fc Yakovlevll2008h and dominates. As a result, r b+s oc (T°°) 2 . At high 
enough temperatures T°° > 6 x 10 8 K the damping due to the bulk viscosity starts to prevail; this results in decreasing of 
Tb+a with growing T°° (the curves in Fig. [4] bend down). At such T°° the neut rons are normal and the proton superfluidity 



is weak or absent. Neglecting the proton superfluidity, one obtains £2 oc (T°°) ijHaensel et al 



l200lT ). hence r b+s oc 1/(T°°) 6 . 
Let us note that the curves for /-modes in Fig. [4] (thick dashed line and thick dots) bend down later than others; for 
them the shear viscosity is the dominant mechanism of damping up to T°° w 2.0 x 10 9 K. This is not surprising, since, 
as it was noted above, for /-modes the deviation from beta-equilibrium is small (8fi°° is reduced by an order of magnitude 
in comparison to p-modes, see Fig. [3]), h ence damping due to the bulk viscosity is suppressed (the relation between Sfi° 
and Tbuik was discussed in detail, e.g., in 
Tb+a oc 1/(T°°) 6 at higher temperatures T 



Gusakov et all 120051 ). As a result, Tb+ S approaches its 'bulk viscosity' asymptote 



> 2.0 x 10 9 K. 



7.3 Frequency spectrum for superfiuid NSs 



First of all let us consider the frequency spectrum for radial oscillations of a superfiuid n eutron star employing the simplified 
model 1 of nucleon superfluidity. For such model this problem was discussed in detail bv lKantor fc Gusakovl (|201ll ). where it 
was solved exactly. Here we compare this exact solution with the approximate calculations obtained in the s = approximation 
(see Sec.[6|. Such a comparison is very useful, since it allows one to make a conclusion about applicability of the approximate 
approach in the case of nonradial oscillations, where the exact solution is not attempted. 

The eigenfrequencies a of radial pulsations (in units of a) versus T$° = T°°/(10 8 K) are shown in Fig.[5ja, b, c). In Fig.[5ja) 
this dependence was obtained assuming that superfiuid and normal modes are completely decoupled (s = approximation). 
The thick solid lines demonstrate the first three normal (nonsuperfluid) radial modes F, 1 H, and 2 H. As one expects, their 
frequencies do not depend on T°° . The dashes are for the first six superfiuid modes 1, . . . , 6, which are the solutions to Eq. 



More precisely, we used Eqs. 1771 and 1781 1 with W s a = V s fl = 0. 
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Figure 5. (color online) Eigenfrcqucncics a (in units of ct) of radial oscillations versus Tg° = T°°/(10 8 K) for model 1 of nucleon 
superfluidity, (a) Approximate spectrum. First three normal modes (F, IH, and 2 H) are shown by the solid lines; first six superfluid 

modes 1, , 6 are shown by dashes, (b) Exact spectrum. Alternate solid and dashed lines show the first six exact modes (I, . . . , VI) of a 

radially oscillating star. No spectrum was plotted in the shaded region, (c) Approximate (dashed lines) and exact (solid lines) spectra. 
At T°° > Tg£ x = 6 X 10 s K all neutrons are normal and the spectrum is that of a nonsuperfluid star. 



These modes, on the contrary, strongly depend on T° 

-,7 



and approach their temperature-independent asymptotes only 



At T°° > T~„ 



at T°° < 5 x 10' K (when the entire NS core is superfluid and Yik does not depend on T 
all neutrons are normal so that superfluid modes do not exist. 

Fig. \Ej[b) demonstrates the results of the exact solution to Eqs. (|83|) - (|101[) obtained by iKantor 



6 x 10 K 



Gusakov (2011) for 



radial oscillations of a superfluid neutron star. The frequencies a of the first six oscillation modes (I,. . .,VI) as functions of 
Tg° are shown by alternate solid and dashed lines. No spectrum is plotted in the gray-shaded area. One can observe that the 
approximate spectrum [Fig.[5ja)] is very similar to the exact spectrum [Fig.[5jb)]. However, there is one important difference: 
instead of crossings of superfluid and normal modes in Fig. EI a) we have avoided crossings of the modes i n Fig. [5Tb). At these 
points the superfluid mode turns into the normal one and vice versa. As it was discussed in details in iGusakov fc Kantoi 
1 201ll ); IKantor &i Gusakov! 1)20111 ). this is not surprising, since in a vicinity of avoided crossings the Einstein equations (|87l) 
and superfluid equation (98]) interact resonantly, so that approximation of completely decoupled superfluid and normal modes 
(s = 0) is inapplicable! 15 . 

For comparison, in Fig. [5jc) we plot both the approximate (dashed lines) and exact (solid lines) spectra. The agreement 
between both spectra is very good: the difference is less than a few per cent. 

Such a close agreement of the exact and approximate results for radial oscillations allows us to analyse the spectrum of 
nonradial oscillations using the same approximation s — 0. The results of this analysis are shown in Fig. [6] for more realistic 
model 2 of nucleon superflu idity (see Se c. [7711 and Fig. [2] ). Superfluid modes shown in this figure have been already studied in 
detail in our recent paper ( Chugunov fc Gusakovll2011l ). Thus, here we discuss them only briefly. 

Fig. [6] contains five panels. Four upper panels present eigenfrequencies a as functions of T%° for normal modes from Table 

1 (thick horizontal lines) and for superfluid modes (dashes) with multipolarities I = 0,1, 2, and 3. For each I there is an infinite 
set of superfluid modes whose eigenfunctions 5/j.i differ by the number of radial nodes n; in the figure we plot the first 25 of them. 
The lower panel demonstrates broadening of the SFL-region with decreasing T£° (SFL-region is shown by hatches). For model 

2 (which we employ here) the redshifted neutron critical temperature T££(r) has a maximum at T^ m ^ « 5.1 x 10 8 K (right 
vertical dotted line). The neutron superfluidity reaches the stellar centre at T°° = T^(0) « 2 x 10 8 K (left vertical dotted line). 
At T°° > T™ max all neutrons are normal, hence only normal modes exist in the star. At T°° < 1^(0) the core is completely 
occupied by the neutron superfluidity. One can see t hat the behaviour of superflu i d modes differs strongly a t T°° > 7^(0) 
and at T°° < T^(0). This feature was discussed in Chugunov fc Gusakov 1 201ll ): Kantor fc Gusakov! l|201ll ). where it was 
demonstrated that (roughly speaking) the frequencies a of superfluid modes scale with Y nn and R B a as a ~ y/Y nn / R s fi, where 
i? s fl is the size of the SFL-region. With the increasing of temperature Y nn decreases, while the size of the SFL-region can 
either decrease [at T°° > T^(0)] or remain constant [at T°° < 7^(0)]. As a result, there is a partial compensation of these 
two tendencies at T°° > 7^(0) (hence, the frequency changes only weakly), while at T°° < 7^(0) the effect of decreasing of 
Ynn is not compensated (hence, a decreases with growing T°°). At T°° < 5 x 10 7 K Yik does not depend on T°°, and, as in 
the case of radial pulsations, the frequencies approach their low-temperature asymptotes. 



Thus, it would not be correct to say that any real oscillation mode of a superfluid star is either purely superfuid or purely normal: 
for some T°° it can show itself as a superfluid, but for other T°° it can behave as a normal mode [see Fig. Ob)] . 
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Figure 6. (color online) Eigenfrequencies <r versus T°° for model 2 of nucleon superfluidity and for multipolarities I = 0, 1, 2, and 3. 
For each I we plot first few normal modes (solid lines) and first 25 superfiuid modes (dashed lines), whose eigenfunctions Sfn differ by 
the number of radial nodes n. At T°° ^ 7^(0) ~ 2 X 10 s K (see the left vertical dotted line), neutron superfluidity occupies the stellar 
centre. The bottom panel demonstrates the variation of the SFL-region (shown by hatching) with T°° . Shaded area in all panels shows 
the region where all neutrons are normal. 
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Figure 7. (color online) Eigenfrequencies a [panels (a, b, c)] and damping times Tb+ S [panels (d, e, f)] of a radially oscillating NS versus 
T°° for model 1 of micleon superfluidity. Panels (a, d): Approximate solution (normal F-mode and first four superfluid modes 1, . . . ,4 are 
shown by solid and dashed lines, respectively); Panels (b, e): Exact solution (first four exact modes I, ... , IV are shown by solid, dashed, 
dot-dashed, and dotted lines, respectively); Panels (c, f): Both approximate (dashed lines) and exact (solid lines) solutions. Panels (a, 
b, c) are the same spectra as those plotted, respectively, in Fig. [5fa, b, c). Normal F-mode is not shown in the shaded region because 
of technical reasons (too many resonances). Dotted lines in panels (d, e, f) show damping times for F-mode calculated using ordinary 
normal-fluid hydrodynamics (see the text for more details). Part of the mode IV (at T°° < 6 X 10 7 K) is shown by dots in the panel (f), 
as described in the text. 



7.4 Damping times for superfluid NSs 

As in the case of eigenfrequencies, we first consider the e-folding times t^ l s = + T~^ eal for radial (Z = 0) pulsations for 
the simplified model 1 of nucleon superfluidity (see Fig. [T]). 

In Fig. [TJa, d) we present the functions cr(T°°) and rb+ s (7 100 ), obtained using the approximate method of Sec. 16.21 The 
frequencies and damping times are plotted for normal F-mode (thick solid line) as well as for the first four superfluid modes 
1, . . . , 4 (dashed lines) [3 In the region shaded in gray the function Tb+a(r°°) for the normal mode was not plotted (there are 
too many merging resonances in this region). The dotted curve in Fig. [TJd, e, f) labeled F n fh ('nfh' is the abbreviation for 
'normal-fluid hydrodynamics') shows the damping time calculated using the ordinary hydrodynamics of nonsuperfluid liquid 
but taking into account the effects of superfluidity on the bulk and shear viscosities. This curve is analogous to the thick solid 
curve in Fig. [4] obtained under the same conditions but for the model 2 of nucleon superfluidity. The vertical dotted line in 
Fig. [7ja, d) indicates a temperature at which frequencies of normal F-mode and the first superfluid mode coincide. 

We present a detailed analysis of Fig. [7jd) in what follows, together with description of the approximate solutions for 
nonradial oscillation modes (Figs. [5] and . 

For comparison, Fig. [T^b, e) demonstrates the results of the exact calculation of frequencies cr(T°°) and damping times 
T h+B (T°°) for the first four (I,. . .,IV) oscillation modes of the superfluid NS [the modes are shown by solid (I), dashed (II), 
dot-dashed (III), and dotted (IV) lines]. 

To see how well the approximate solution [Fig.[7Ja, d)] agrees with the exact one [Fig.[7Jb, e)], both solutions are presented 
in Fig. [7Jc, f). Dashes correspond to approximate solution, solid lines - to exact solution. A portion of the mode IV in Fig. 
[TJf) is shown by thick dots because the corresponding approximate solution (the mode 1 H) is not plotted. One sees that the 
agreement between the approximate and exact solutions is reasonable everywhere (average error does not exceed 10 — 25%) 
except for the resonances (see below) and an interval of temperatures T°° < 3 x 10 7 K where the mode III of exact solution 



16 Notice that, in Figs. [3a, b, c) we present, in logarithmic scale, parts of the spectra, which were already plotted in linear scale in Figs. 
(5fa, b, c), respectively. 
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Figure 8. (color online) Damping times Tb-|_ a versus T°° for various oscillation modes for model 2 of nuclcon superfluidity. On each 
panel we plot one normal mode (shown by solid line; its multipolarity and name is indicated) and first f 5 superfluid modes (dashed 
lines). Dotted lines show r b _|_ s (T°°) for normal modes calculated using normal-fluid hydrodynamics [see the text for more details]. In the 
shaded area all neutrons are normal and superfluid modes do not exist. 



deviates from the second superfluid mode of approximate solution. To explain this deviation let us note that, as follows from 
Fig. [5{a), at such T°° the frequency of the normal mode 1H practically coincides with that of the second superfluid mode. 
In that case Eqs. (|87[) and (|98l) interact resonantly, so that the approximation of independent superfluid and normal modes 
is poor even though parameter s is small 

Let us now consider the nonradial oscillations. Fig. [5] presents an approximate solution for the function Tb+ S (r°°), which 
is obtained for a realistic nucleon superfluidity model 2. By dashes we show superfluid modes, solid lines correspond to normal 
modes. Each panel in the figure is plotted for one normal mode (its name and multipolarity / are indicated) and for the first 
15 superfluid modes with the same I. By dots, as in Fig. [7^d, e, f), we plot Tb+ S for a corresponding normal modes calculated 
using the ordinary normal-fluid hydrodynamics. In the shaded region superfluid modes were not plotted because all neutrons 
are normal there and the star oscillates as a nonsuperfluid. 

In more detail damping times are demonstrated for quadrupole (/ = 2) oscillation modes in Fig. [9] In particular, the 
normal pi-mode is shown there by solid lines. In the three lower panels we plot the dependence Tb+s(T°°) in an increasingly 
larger scale. In the three upper panels we plot, in the same scale, the oscillation frequencies a(T°°) (the corresponding spectrum 
was already presented in Fig. [6] in linear scale). Left lower panel of Fig. [9] coincides with Fig. [8{e). 

Let us discuss the main conclusions that can be drawn from the analysis of Figs. [7{d), [8] and [9] 

1. For any normal mode the dependence Tb+ a (T°°) (solid lines in these figures), has a set of resonance features (spikes) 
concentrated (for radial and p- modes) to the critical temperature 7^(0) at which neutron superfluidity in the core centre 
dies out. For model 1 1^(0) = T~ max = 6 x 10 8 K (see Fig. [TJ) , for model 2 T c ^(0) w 2 x 10* K (see Fig.©. The resonances 
appear when frequency of the normal mode approaches the frequency of one of the superfluid modes. For instance, solid line 
in Fig. [7J a ) crosses superfluid modes four times [in Fig. [TJa, d) the temperature T°° of the first crossing is shown by the 



17 For a quantitatively correct description of the function ib+ s (T°°) in that case it is, in principle, straightforward to develop a pertur- 
bation theory similar to the degenerate perturbation theory of quantum mechanics (see below the discussion of resonances in Figs. [7] El 
and 
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Figure 9. (color online) Eigenfrequcncies <r (upper panels) and damping times Tb-|_ s (lower panels) versus T°° for quadrupole (I = 2) 
oscillation modes in an increasingly larger scale. The normal pi-mode is shown by solid lines. Left lower panel coincides with Fig. |8je). 
Other lower panels are zoomed in versions of Fig. [8[e) . Notations are the same as in Figs. [6] and [8] 

vertical dotted line and equals T°° w 10 s K]. Correspondingly, four resonances appear in Fig. [7Jd). A similar situation can 
be observed in Figs. [8] and [9] Near resonances Tb+ S for normal mode rapidly decreases by 1-2 orders of magnitude (see item 
2 below) and, in the resonance point, it becomes strictly equal to Tb+ S for the corresponding superfluid mode. 

Such behavior of the approximate solution Tb+ S (r°°) for normal modes in the vicinity of resonances can be easily un- 
derstood. In resonance points, in which the frequencies of superfluid and normal modes coincide, Eq. (|98l) has a nontrivial 
solution even in the absence of the source S/j, n oimi- For it to be satisfied with the source, the oscillation amplitude Sfj, must 
be infinitely large. In other words, in resonance points all the energy must be contained in superfluid degrees of freedom (in 
particular, near resonances W s a S> W4> and V B n 3> 14,). Formally, this means that in the resonance point the damping time 
rb+s should be exactly the same as for the superfluid mode. 

Another important point that is worth noting is that, as follows from Fig. [7jf), the approximate solution for the normal 
radial F-mode describes qualitatively well the exact solution near resonances (the latter is shown by solid lines). We expect 
that the same is also true for nonradial modes for which the exact solution was not attempted. At first glance such an 
agreement between the approximate and exact solutions seems surprising because the approximation s = should not work 
in the vicinity of resonances, where the frequencies of superfluid and normal modes are close to each other. Nevertheless, 
one verifies that this approximation is still suitable for a qualitatively correct description of the function Tb+ S (T°°) if one 
bears in mind that: (i) close to any resonance the exact solution is a linear superposition of independent solutions describing 
(intersecting) superfluid and normal modes and (ii) Tb+ S for the superfluid mode is much less than for the normal mode. 

Items (i) and (ii) mean that, in the exact solution, the main contribution to Tb+ S comes from the superfluid mode 
(while the contribution from the normal mode is small). This leads us to conclusion that the superfluid modes are the main 
sources of viscous dissipation in the vicinity of resonance points. The same conclusion was already drawn above using the 
approximate method of Sec. 16.21 This explains why the approximate method gives qualitatively correct results for Tb+ S {T°°) 
near resonances. 

In order to avoid confusion let us emphasize that the function Tb+ S (T°°) contains resonance features (spikes) for normal 
modes only in the approximate solution [see Figs. [7Jd),[8] and [9]. In the exact solution any normal oscillation mode turns 
into a superfluid one near resonance (and vice versa). This leads to an abrupt decreasing (increasing) of 7~b+ s and formation 
of a 'step-like' structure rather than spike [see Fig.|7[e)]. 

2. It was already mentioned above that, as follows from Figs. [TJd), [8] and [9] normal modes (far from resonances) damp 
out by 1-2 orders of magnitude slower than those superfluid modes with which they can have equal frequencies (i.e. intersect 
in the a - T°° plane). 

What is the reason for such a fast damping of superfluid modes? To be more concrete, below we consider a low-temperature 
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case, T°° < 3 x 10 7 K. There are three main factors: (i) For superfluid modes eigenfunctions W B n and V B a have a maximum 
in the central regions of a star where the shear viscosity is maximal. On the contrary, for normal modes the maximum of 
eigenfunctions Wb and Vb lies closer to the NS surface, where the shear viscosity coefficient can be substantially (5 and more 
times) smaller. As a consequence, 2H s hcar for superfluid modes turns out to be greater (and hence r s h oa r smaller) than for 
normal modes, (ii) The energy of superfluid modes is given by Eq. (|73[) and depends on the quantity y [see Eq. (I68p for 
the definition of y]. At low T°° the parameter y is small, y ~ n p /n n ~ 0.04 -j- 0.09, which also results in decreasing of the 
characteristic damping times for superfluid modes 0- (Hi) This factor is particularly important for radial oscillations (I — 0) 
and is related to a coefficient a± in the expression (|78[) for the damping time T s hcar due to shear viscosity. This coefficient is 
given by Eq. (|8f | which is a sum of four terms. It turns out that for the normal radial modes the first term is well compensated 
by the third term H2, while the other therms vanish. For the superfluid modes such compensation does not occur because for 
them H 2 = 0. 

3. At low enough T°° the damping times for normal radial and p- modes can be several times larger or smaller that Tb+ S , 
calculated employing ordinary hydrodynamics of nonsuperfluid liquid but accounting for the effects of superfluidity on the 
bulk and shear viscosity coefficients (dotted lines in Figs. [7jd), |8ja, d, e, f), and [9}. Let us inspect, for example, Fig. [8ja). 
One sees that at T°° < 10 8 K Tb+ S , calculated in the frame of nonsuperfluid hydrodynamics, is approximately 4 times larger 
than rb+ s determined self-consistently. This difference arises because to plot the dotted curve we used the formulas of Sec. |4] 
in which W s a = V a n — 0. As T°° grows, however, the difference in two ways of calculating Tb+ S rapidly decreases because the 
SFL-region becomes smaller and hence its contribution to Tb+ S becomes less and less pronounced. 

4. Unlike the radial and p-modes, the agreement between dotted and solid lines for normal /-modes is very good [see Fig. 
[8jb, c)], which means that for these modes use of the nonsuperfluid hydrodynamics (far from resonances) is well justified. 
The reason for such a good agreement of damping times is related to a relatively weak compression-decompression of matter 
in the course of the /-type oscillations. As a consequence, for the normal /-modes the source <5/x norm ; in Eq. (|98[1 is small, so 
that far from the resonances S^j, ~ and the superfluid degrees of freedom are almost not excited [W s a ~ Kfl ~ 0, see Eqs. 
(|48|l . (|52[) . and (|90[) ]. This result confirms, e xtends and, we think, provid es a deeper understanding, o f the results previously 
obtained in a Newtonian framework by, e.g., Lindblom fc Mendell (1994) and Andersson et al. 1 20091 ). 



5. At T°° — > T^ max = 6 x 10 8 K one can observe the rapid increasing of Tb+ S for superfluid modes in Fig. [7] It is bounded 
from above by Tbuik and is related with the ten dency of T snca r to grow to infinity in this limit. Such a behaviour of T sn0 ar was 
discussed in detail in lKantor fc Gusakovl (|201 lh and is specific for model 1 of nucleon superfluidity. 



8 SUMMARY 



In this paper we, for the first time, self-consistently analyze the effects of nucleon superfluidity on damping of oscillations of 
nonrotating general relativtstic NSs. Our main results are summarized below. 

1. The analytic formulas are derived for the oscillation energy -E mcc h (|7ip and for the characteristic damping times Tbuik 
(|77|l and r s h oa r (|78|l due to the bulk and shear viscosities. These expressions ar e valid for oscillations of a r bitrar y multipolarity 
I. The expression (|71[1 for -E mccn is the generalization of the formula (26) of iThorne fc Campolattard |l967), written for a 
nonsu perfluid NS. The ex pressions (|77|l and (|78p are the generalizations, to the case of sup erfluidity, of the for mulas (5) and 
(6) in lCutler et al.1 1 1990l ). Notice that the damping times, calculated using the formulas of Cutler et al. ( 1990h appear to be 
2 times smaller than our Tbuik and Tshear, calculated from Eqs. (|77[) and (|78[l under assumption that superfluid degrees of 
freedom are suppressed (i.e., W B a = V s n = 0). 

2. An approximate method is developed in detail and applied, which allows one to easily determine the eigenfrequencies 
and eigenfunctions of an oscillating superfluid NS, provided that they are known for a normal (nonsuperfluid) star of the 
same mass (see Sec. I6.2JI . The method is based on the appr oximate decoupling of equ a tions describing superfluid an d normal 
oscillation modes and exploits the ideas first formulated in Gusakov fc Kanton ( 2011 ): Chueunov fc Gusakov 1 201lh . 

3. Using radial oscillations as an example, and adopting the simplified model 1 of nucleon superfluidity (Fig. (TJ, we 
demonstrate that this method leads to oscillation frequencies and characteristic damping times that agree well with the 
results of exact calculation. 

4. The approximate method of Sec. l6.2l is applied to study nonradial oscillations of a superfluid NS assuming the realistic 
model 2 of nucleon superfluidity (Fig. [2}. A number of normal and superfluid oscillation modes with multipolarities I = 0, ... ,3 
are considered. In particular, the following normal modes are analyzed: F-mode for I = 0, pi-mode for I — 1, f- and pi-modes 
for I — 2 and 3. 

It is demonstrated that: 

(i) As a rule, for any given normal mode (whose frequency a coincides with the corresponding frequency of a nonsuperfluid 



8 To get an estimate for y we made use of the sum rule u n Y nn + fi p Y np 
and neglected the small matrix element Y np in comparison with In- 



valid at T°° = (Gu sakov. Kantor & Hacnscl 2009a), 
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NS and does not depend on the internal redshifted stellar temperature T°°) the viscous damping time Tb+ S = ( 7 "buik + r shear^ 
is one order of magnitude greater than t^ +b for those superfluid modes that can intersect the normal mode in the a — T°° 
plane. This effect is non-local (occurs only after integration over the NS volume) and is determined by a number of factors 
(see item 2 of Sec. I7.4[>. 

(ii) The function Tb+ S (T°°) for any normal mode contains resonance features. In resonance points the frequency a of a 
normal mode coincides with that of some of the superfluid modes (their a depend on T°°). When passing a resonance (e.g., 
with growing T°°), Tb+ S initially rapidly decreases (by 1-2 orders of magnitude) until it reaches the value of Tb+ S for this 
superfluid mode and, after that, it increases again (see Figs. [7jd), [8j and[9|. 

(Hi) Resonance features (spikes) appear only in the approximate treatment of Sec. [6] in which the normal and superfluid 
modes intersect at resonance points [see, e.g., Fig. Ufa)]. In the exact solution instead of crossings one has avoided crossings 
of modes [Fig. [Sfb)]. Near avoided crossings any real mode changes its behaviour from normal- like to superfluid- like (and vice 
versa). As a result, instead of spikes one has a very rapid step- like decreasing (increasing) of Tb+ S [cf. Figs. [TJd) and[7Je)]. 

(iv) Sufficiently far from the resonances 7~b+ s for normal radial and p-modes, determined self-consistently employing the 
hydrodynamics of a superfluid liquid, can differ several fold from 7b+ s , calculated using the ordinary normal-fluid hydrody- 
namics (but accounting for the effects of superfluidity on the shear and bulk viscosities). The latter approximation is often 
adopted in the literature devoted to oscillations of NSs. 

(v) In contrast to radial and p-modes, for /-modes far from the resonances, use of the ordinary hydrodynamics of 
nonsuperfluid liquid for calculation of Tb+ S is well justified. The reason is that for /-type oscillations the imbalance S/j, of 
chemical potentials is relatively small (matter does not compress significantly during oscillations). Thus, superfluid degrees of 
freedom are almost not excited (see Sees. 16.21 and [774)) . 

(vi) Since for /-modes far from the resonances S/j, is small (that is, deviation from the beta-equilibrium is weak), bulk 
viscous damping of /-modes is suppressed in comparison to p-modes. 

Though here we only considered oscillations of superfluid nonrotating NSs, we expect that the main conclusions of this 
work will also remain (mostly) unchanged for rotating NSs. Our results indicate that dissipative evolution of oscillating NSs 
may follow quite different scenarios than those usually considered in the literature. This is especially true if one is interested 
in the combined analysis of damping of oscillations and thermal evolution of a NS or in the analysis of instability windows, 
that is the values of T°° and rotation fr equency at which a star becomes unstable wit h respect to the emission of gravitational 
waves (e.g., the r-mode instability, see Andersson 19981 ; Friedman &: Morsink 1998T ). These issues are extremely interesting 
and important, but we left them beyond the scope of the present paper and will address the related topics in our subsequent 
publication. 
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APPENDIX A: BOUNDARY CONDITIONS TO EQUATION (98) 

Equation (|98[) should be solved in the region of a NS core where neutrons are superfluid (SFL-region). If the NS centre is 
occupied by the neutron superfluidity, then for regularity of the solution at r — > it is necessary that 

Sfii oc r . (Al) 

The conditions at the boundary of the SFL-region follow from the requirement of the absence of particle transfer (baryons 
and electrons) through the interface. One obtains from the definitions ©-([7]) 

X ± = 0, (A2) 

where X± is the component of the vector X J perpendicular to the interface. To rewrite Eq. (|A2[) in terms of 5(ii(r), it is 
necessary to consider two possibilities: 

(i) The boundary (one of the boundaries) between the SFL-region and nonsuperfluid matter lies inside the core and 
is defined by the condition T = T cn (Rb) [Rb is the radial coordinate of the boundary]. Then at the boundary Y na (Rb) ~ 
Y n p(Rb) = and from Eqs. ([90]) and (J98J one has 

e x ' v/2 w 2 

8^1 = TT^a _ Sfimmal). (A3) 
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(ii) Outer boundary of the SFL-region coincides with the crust-core interface (i?b 
[that is Ynn(-Rcc) and Y np (-R C c) are non-zero] and from Eq. (|90[) it follows that 

<5/4(i*cc) =0. 

The conditions (|A1[) — <|A4[) are necessary and sufficient for solving Eq. (|98[) . 
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